home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / Tools / Dev / SpeakFreely_Src / gsm / src / lpc.c < prev    next >
C/C++ Source or Header  |  2000-05-27  |  7KB  |  344 lines
  1. /*
  2.  * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
  3.  * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
  4.  * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
  5.  */
  6.  
  7. /* $Header: /home/kbs/jutta/src/gsm/gsm-1.0/src/RCS/lpc.c,v 1.1 1992/10/28 00:15:50 jutta Exp $ */
  8.  
  9. #include <stdio.h>
  10. #include <assert.h>
  11.  
  12. #include "private.h"
  13.  
  14. #include "gsm.h"
  15. #include "proto.h"
  16.  
  17. #undef    P
  18.  
  19. /*
  20.  *  4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
  21.  */
  22.  
  23. /* 4.2.4 */
  24.  
  25.  
  26. static void Autocorrelation P2((s, L_ACF),
  27.     word     * s,        /* [0..159]    IN/OUT  */
  28.      longword * L_ACF)    /* [0..8]    OUT     */
  29. /*
  30.  *  The goal is to compute the array L_ACF[k].  The signal s[i] must
  31.  *  be scaled in order to avoid an overflow situation.
  32.  */
  33. {
  34.     register int    k, i;
  35.  
  36.     word        temp, smax, scalauto;
  37.     longword    L_temp;
  38.  
  39. #ifdef    USE_FLOAT_MUL
  40.     float        float_s[160];
  41.     longword    L_temp2;
  42. #else
  43.     longword    L_temp2;
  44. #endif
  45.  
  46.     /*  Dynamic scaling of the array  s[0..159]
  47.      */
  48.  
  49.     /*  Search for the maximum.
  50.      */
  51.     smax = 0;
  52.     for (k = 0; k <= 159; k++) {
  53.         temp = GSM_ABS( s[k] );
  54.         if (temp > smax) smax = temp;
  55.     }
  56.  
  57.     /*  Computation of the scaling factor.
  58.      */
  59.     if (smax == 0) scalauto = 0;
  60.     else {
  61.         assert(smax > 0);
  62.         scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
  63.     }
  64.  
  65.     /*  Scaling of the array s[0...159]
  66.      */
  67.  
  68.     if (scalauto > 0) {
  69.  
  70. # ifdef USE_FLOAT_MUL
  71. #   define SCALE(n)    \
  72.     case n: for (k = 0; k <= 159; k++) \
  73.             float_s[k] = (float)    \
  74.                 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
  75.         break;
  76. # else 
  77. #   define SCALE(n)    \
  78.     case n: for (k = 0; k <= 159; k++) \
  79.             s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
  80.         break;
  81. # endif /* USE_FLOAT_MUL */
  82.  
  83.         switch (scalauto) {
  84.         SCALE(1)
  85.         SCALE(2)
  86.         SCALE(3)
  87.         SCALE(4)
  88.         }
  89. # undef    SCALE
  90.     }
  91. # ifdef    USE_FLOAT_MUL
  92.     else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
  93. # endif
  94.  
  95.     /*  Compute the L_ACF[..].
  96.      */
  97.     {
  98. # ifdef    USE_FLOAT_MUL
  99.         register float * sp = float_s;
  100.         register float   sl = *sp;
  101. # else
  102.         word  * sp = s;
  103.         word    sl = *sp;
  104. # endif
  105.  
  106. #    define STEP(k)     L_ACF[k] += (longword)(sl * sp[ -(k) ]);
  107. #    define NEXTI     sl = *++sp
  108.  
  109.  
  110.     for (k = 9; k--; L_ACF[k] = 0) ;
  111.  
  112.     STEP (0);
  113.     NEXTI;
  114.     STEP(0); STEP(1);
  115.     NEXTI;
  116.     STEP(0); STEP(1); STEP(2);
  117.     NEXTI;
  118.     STEP(0); STEP(1); STEP(2); STEP(3);
  119.     NEXTI;
  120.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
  121.     NEXTI;
  122.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
  123.     NEXTI;
  124.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
  125.     NEXTI;
  126.     STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
  127.  
  128.     for (i = 8; i <= 159; i++) {
  129.  
  130.         NEXTI;
  131.  
  132.         STEP(0);
  133.         STEP(1); STEP(2); STEP(3); STEP(4);
  134.         STEP(5); STEP(6); STEP(7); STEP(8);
  135.     }
  136.  
  137.     for (k = 9; k--; L_ACF[k] <<= 1) ; 
  138.  
  139.     }
  140.     /*   Rescaling of the array s[0..159]
  141.      */
  142.     if (scalauto > 0) {
  143.         assert(scalauto <= 4); 
  144.         for (k = 160; k--; *s++ <<= scalauto) ;
  145.     }
  146. }
  147.  
  148. #if defined(USE_FLOAT_MUL) && defined(FAST)
  149.  
  150. static void Fast_Autocorrelation P2((s, L_ACF),
  151.     word * s,        /* [0..159]    IN/OUT  */
  152.      longword * L_ACF)    /* [0..8]    OUT     */
  153. {
  154.     register int    k, i;
  155.     float f_L_ACF[9];
  156.     float scale;
  157.  
  158.     float          s_f[160];
  159.     register float *sf = s_f;
  160.  
  161.     for (i = 0; i < 160; ++i) sf[i] = s[i];
  162.     for (k = 0; k <= 8; k++) {
  163.         register float L_temp2 = 0;
  164.         register float *sfl = sf - k;
  165.         for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
  166.         f_L_ACF[k] = L_temp2;
  167.     }
  168.     scale = MAX_LONGWORD / f_L_ACF[0];
  169.  
  170.     for (k = 0; k <= 8; k++) {
  171.         L_ACF[k] = f_L_ACF[k] * scale;
  172.     }
  173. }
  174. #endif    /* defined (USE_FLOAT_MUL) && defined (FAST) */
  175.  
  176. /* 4.2.5 */
  177.  
  178. static void Reflection_coefficients P2( (L_ACF, r),
  179.     longword    * L_ACF,        /* 0...8    IN    */
  180.     register word    * r            /* 0...7    OUT     */
  181. )
  182. {
  183.     register int    i, m, n;
  184.     register word    temp;
  185.     register longword ltmp;
  186.     word        ACF[9];    /* 0..8 */
  187.     word        P[  9];    /* 0..8 */
  188.     word        K[  9]; /* 2..8 */
  189.  
  190.     /*  Schur recursion with 16 bits arithmetic.
  191.      */
  192.  
  193.     if (L_ACF[0] == 0) {    /* everything is the same. */
  194.         for (i = 8; i--; *r++ = 0) ;
  195.         return;
  196.     }
  197.  
  198.     assert( L_ACF[0] != 0 );
  199.     temp = gsm_norm( L_ACF[0] );
  200.  
  201.     assert(temp >= 0 && temp < 32);
  202.  
  203.     /* ? overflow ? */
  204.     for (i = 0; i <= 8; i++) ACF[i] = SASR( L_ACF[i] << temp, 16 );
  205.  
  206.     /*   Initialize array P[..] and K[..] for the recursion.
  207.      */
  208.  
  209.     for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
  210.     for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
  211.  
  212.     /*   Compute reflection coefficients
  213.      */
  214.     for (n = 1; n <= 8; n++, r++) {
  215.  
  216.         temp = P[1];
  217.         temp = GSM_ABS(temp);
  218.         if (P[0] < temp) {
  219.             for (i = n; i <= 8; i++) *r++ = 0;
  220.             return;
  221.         }
  222.  
  223.         *r = gsm_div( temp, P[0] );
  224.  
  225.         assert(*r >= 0);
  226.         if (P[1] > 0) *r = -*r;        /* r[n] = sub(0, r[n]) */
  227.         assert (*r != MIN_WORD);
  228.         if (n == 8) return; 
  229.  
  230.         /*  Schur recursion
  231.          */
  232.         temp = GSM_MULT_R( P[1], *r );
  233.         P[0] = GSM_ADD( P[0], temp );
  234.  
  235.         for (m = 1; m <= 8 - n; m++) {
  236.             temp     = GSM_MULT_R( K[ m   ],    *r );
  237.             P[m]     = GSM_ADD(    P[ m+1 ],  temp );
  238.  
  239.             temp     = GSM_MULT_R( P[ m+1 ],    *r );
  240.             K[m]     = GSM_ADD(    K[ m   ],  temp );
  241.         }
  242.     }
  243. }
  244.  
  245. /* 4.2.6 */
  246.  
  247. static void Transformation_to_Log_Area_Ratios P1((r),
  248.     register word    * r             /* 0..7       IN/OUT */
  249. )
  250. /*
  251.  *  The following scaling for r[..] and LAR[..] has been used:
  252.  *
  253.  *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
  254.  *  LAR[..] = integer( real_LAR[..] * 16384 );
  255.  *  with -1.625 <= real_LAR <= 1.625
  256.  */
  257. {
  258.     register word    temp;
  259.     register int    i;
  260.  
  261.  
  262.     /* Computation of the LAR[0..7] from the r[0..7]
  263.      */
  264.     for (i = 1; i <= 8; i++, r++) {
  265.  
  266.         temp = *r;
  267.         temp = GSM_ABS(temp);
  268.         assert(temp >= 0);
  269.  
  270.         if (temp < 22118) {
  271.             temp >>= 1;
  272.         } else if (temp < 31130) {
  273.             assert( temp >= 11059 );
  274.             temp -= 11059;
  275.         } else {
  276.             assert( temp >= 26112 );
  277.             temp -= 26112;
  278.             temp <<= 2;
  279.         }
  280.  
  281.         *r = *r < 0 ? -temp : temp;
  282.         assert( *r != MIN_WORD );
  283.     }
  284. }
  285.  
  286. /* 4.2.7 */
  287.  
  288. static void Quantization_and_coding P1((LAR),
  289.     register word * LAR        /* [0..7]    IN/OUT    */
  290. )
  291. {
  292.     register word    temp;
  293.     longword    ltmp;
  294.     ulongword    utmp;
  295.  
  296.  
  297.     /*  This procedure needs four tables; the following equations
  298.      *  give the optimum scaling for the constants:
  299.      *  
  300.      *  A[0..7] = integer( real_A[0..7] * 1024 )
  301.      *  B[0..7] = integer( real_B[0..7] *  512 )
  302.      *  MAC[0..7] = maximum of the LARc[0..7]
  303.      *  MIC[0..7] = minimum of the LARc[0..7]
  304.      */
  305.  
  306. #    undef STEP
  307. #    define    STEP( A, B, MAC, MIC )        \
  308.         temp = GSM_MULT( A,   *LAR );    \
  309.         temp = GSM_ADD(  temp,   B );    \
  310.         temp = GSM_ADD(  temp, 256 );    \
  311.         temp = SASR(     temp,   9 );    \
  312.         *LAR  =  temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
  313.         LAR++;
  314.  
  315.     STEP(  20480,     0,  31, -32 );
  316.     STEP(  20480,     0,  31, -32 );
  317.     STEP(  20480,  2048,  15, -16 );
  318.     STEP(  20480, -2560,  15, -16 );
  319.  
  320.     STEP(  13964,    94,   7,  -8 );
  321.     STEP(  15360, -1792,   7,  -8 );
  322.     STEP(   8534,  -341,   3,  -4 );
  323.     STEP(   9036, -1144,   3,  -4 );
  324.  
  325. #    undef    STEP
  326. }
  327.  
  328. void Gsm_LPC_Analysis P3((S, s,LARc),
  329.     struct gsm_state *S,
  330.     word          * s,        /* 0..159 signals    IN/OUT    */
  331.         word          * LARc)    /* 0..7   LARc's    OUT    */
  332. {
  333.     longword    L_ACF[9];
  334.  
  335. #if defined(USE_FLOAT_MUL) && defined(FAST)
  336.     if (S->fast) Fast_Autocorrelation (s,      L_ACF );
  337.     else
  338. #endif
  339.     Autocorrelation              (s,      L_ACF    );
  340.     Reflection_coefficients          (L_ACF, LARc    );
  341.     Transformation_to_Log_Area_Ratios (LARc);
  342.     Quantization_and_coding          (LARc);
  343. }
  344.